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Abstract 

We present a new technique for the numerical simulation of axisymmetric 
systems. This technique avoids the coordinate singularities which often arise 
when cylindrical or polar-spherical coordinate finite difference grids are used, 
particularly in simulating tensor partial differential equations like those of 
3 + 1 numerical relativity. For a system axisymmetric about the z axis, the 
basic idea is to use a 3-dimensional Cartesian (x, y, z) coordinate grid which 
covers (say) the y = plane, but is only one finite-difference-molecule-width 
thick in the y direction. The field variables in the central y = grid plane can 
be updated using normal (x, y, ^-coordinate finite differencing, while those 
in the y 7^ grid planes can be computed from those in the central plane by 
using the axisymmetry assumption and interpolation. We demonstrate the 
effectiveness of the approach on a set of fully nonlinear test computations in 
3 + 1 numerical general relativity, involving both black holes and collapsing 
gravitational waves. 
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I. INTRODUCTION 



Finite difference numerical simulations of axisymmetric systems are most often, and most 
naturally, carried out in coordinate systems adapted to the symmetry of the underlying 
problem, e.g. polar spherical (r, 9, <fi) or cylindrical (p, z, <fi) coordinates. However, the use of 
such coordinate systems brings with it delicate problems in finite differencing near the z axis, 
particularly when tensor time-evolution partial differential equations (PDEs) are considered. 
Depending on the problem, it is often very difficult to obtain fully stable numerical evolutions 
near the axis, and for some problems it is difficult to even accurately discretize the equations 
there. Here we consider the problem for the Einstein equations of general relativity in 
axisymmetry, although our approach should be useful for other systems of PDEs (e.g. the 
Navier-Stokes equations), and also for the case of spherical symmetry. 

There are several different types of z axis difficulties, which depending on the physi- 
cal system may occur singly or in combination. The simplest problem is that physically 
nonsingular terms in the equations may have indeterminate 0/0 forms along the z axis. 
Fortunately, such terms are generally easy to regularize by applying L'Hopital's rule. For 
example, in polar spherical coordinates (r,9,(f>), the flat-space Laplacian operator includes 
the term 



1 —d e ( S m9d e ). (1) 



r 2 sin 9 

Assuming all fields to be smooth on the z axis and applying L'Hopital's rule, the 9 — > 
limit of this term is easily seen to be 

y 2 d ee . (2) 

However, even here there may be a problem with finite differencing: Although the 9 — ► 
limit of the term ([!]) is precisely (Q), when we finite difference these terms the limiting 
relationship generally only holds in the limit where the grid spacing A9 — > 0. For any given 
(nonzero) A9, the numerically computed values of ([]]) will have a 9 —>■ limit which will in 
general differ somewhat from the numerical values of (fj). This difference may give rise to 
finite differencing instabilities near the axis. 

A more serious z axis problem is that of oo — oo cancellations, where a physically nonsin- 
gular quantity is computed as the sum of many terms, two or more of which are individually 
singular on the z axis. For example, again using polar spherical coordinates, (r,9,(f)), in 
the so-called "3 + 1" Einstein equations ( see, for example, |2j|| for general reviews), 
the 3-Ricci tensor component Rgg is generically 0(1) on the z axis, but it is computed as 
the sum of, among (many) other terms, —^g^dggg^, which is generically 0(l/9 2 ) near the 
z axis. Although not impossible, regularizing terms of this nature is very difficult, requiring 
a detailed analysis of the generic behavior of the entire system of PDEs under consideration 
(for our example, the coupled 3 + 1 evolution of Einstein's equations, with constraint and 
coordinate equations) near the z axis. 

Collectively these problems can be very severe, crippling attempts to evolve systems with 
such symmetries (see, e.g. [§-0]). Many solutions to the problems brought on by special co- 
ordinate systems have been attempted, including various regularization procedures fg-11 
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Taylor expansions ||10|| , the use of nonsingular-basis tensor components 0, spectral meth- 
ods [12 1, and special coordinate conditions used to eliminate troublesome terms [|,13,5|- 



However, these methods tend to be complicated and not particularly robust. 

On the other hand, 3D Cartesian coordinates contain no coordinate pathologies, and all 
terms in the equations governing the evolution of functions are typically completely regular, 
even at special points such as the origin or along a symmetry axis |14]|. Normally, however, 
if one treats a spherical or axisymmetric system in 3D Cartesian coordinates one loses the 
ability to ignore irrelevant dimensions. For example, a spherically symmetric system should 
reduce to a ID problem, depending only on the radius r, while an axisymmetric system 
reduces to a 2D problem, independent of the azimuthal coordinate <fi. In full 3D Cartesian 
coordinates, not only is the true symmetry of the problem potentially disturbed by the 
finite differences taken in a Cartesian coordinate system, but dimensional reductions do 
not occur, and the memory requirements are much larger, scaling as N 3 , rather than as 
iV 2 or N. For reasonable grid sizes of order hundreds of zones, these factors can become 
astronomical. In 3D numerical relativity, the problem is exacerbated by the need to carry 
a large number of 3D grid functions (typically over 100) in memory at all times. Some 
savings have been realized in cases where Cartesian coordinates are used for intrinsically 
axisymmetric or spherical problems by evolving only a single octant, reducing memory and 
computational requirements by a factor of eight, but this does not change the overall scaling 
with N. 

Furthermore, while many problems are now treatable in 3D, where Cartesian coordinates 
are generally favored, testbeds for developing algorithms for the full 3D case are often carried 
out in lower dimensional cases, where the problems are simpler and where limiting solutions 
are known. However, working in special coordinate systems for lower dimensional test 
problems introduces difficulties specific to the coordinate system, and frequently techniques 
developed in special coordinate system do not carry over to the 3D Cartesian code. What 
is needed is a lower dimensional testbed that retains the same essential features present in 
the generic 3D case, at both the physical and computational level. 

In this paper we describe a scheme which borrows from the singularity-free nature of 
full 3D Cartesian coordinates, and allows the treatment of axisymmetric and spherically 
symmetric problems without the memory constraints of full 3D Cartesian coordinates. Take 
the case of axisymmetry: The essential trick is to realize that in 3D Cartesian coordinates, 
an axisymmetric system can be computed in, say, the x-z (y = 0) plane alone. The x-z 
system can be rotated about the z axis to determine the solution at any (x, y, z) point 
at a given instant of time. However, a 3D evolution also requires spatial derivatives in the 
y direction, which (for non-scalar field variables) do not necessarily all vanish, even in the x-z 
(y = 0) plane. Because of the axisymmetry assumption, the solution in the x-z (y = 0) plane 
can be rotated according to tensor transformation laws to any y value, so the y derivatives 
of all quantities can be determined in the x-z plane solely from information in this same 
plane. Hence, a full 3D evolution in Cartesian coordinates can be carried out, using only 
information from a single 2D plane of data. There are different ways to achieve this, as we 
show in the sections below, but the key point is that an axisymmetric or spherical system 
can be evolved as if it were in 3D Cartesian coordinates , and hence without coordinate 
singularities and the instabilities they can induce, and without going to the expense of a full 
3D computation. 
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In section |I| we describe the technique in detail for the case of axisymmetry. We show 
the effectiveness of the method in applications to dynamic wave and black hole spacetimes in 
section HI, and summarize the work in section [TV]. The method has been fully implemented 
under the name "Cartoon" (chosen because of its resemblance to the way low-budget tele- 
vision cartoons animate a nominally 3-dimensional world in "2~" dimensions, and also as 
a shorthand for cartesian two-dimensional) in the Cactus code for numerical relativity and 
astrophysics ]ToVP^l . 



II. THE TECHNIQUE 

In 3D systems with discrete symmetries it has been common practice to evolve just a part 
of the domain, and use symmetries to provide boundary conditions. For example, spherical, 
axisymmetric, and even 3D systems have been evolved in numerical relativity in a single 
octant, if the symmetries of the system allow it. A spherical system is reflection symmetric 
about all coordinate planes, and hence it could be evolved in Cartesian coordinates in the 
octant defined by, say, x, y, z all non-negative. Any planar reflection symmetry can be used 
to provide boundary conditions at the plane of reflection, i.e. from the symmetry one can 
infer that certain functions are either symmetric or antisymmetric across the coordinate 
planes. This has been used in many simulations in numerical relativity |I7| , |T9| -|26 . 



In what follows we discuss how in the case of axisymmetry one can use the continu- 
ous rotational symmetry to provide boundary conditions for a thin 3D slab. The same 
construction is expected to be applicable to spherically symmetric systems. 

Consider a rectangular three dimensional volume in R 3 , with Cartesian coordinates 
(x,y,z). We assume a uniformly spaced Cartesian grid with spacings Ax, Ay, and Az, 
and we take the finite difference molecules to have radii v in the x direction and w in the 



y direction. (For reasons discussed in section fl~B[ we will typically have v > w.) We take 



the z axis to be the axis of symmetry, as shown in Fig. [T]. In the y direction, the grid con- 
tains 2w + 1 points centered about y = 0. (For example, for standard second order centered 
difference w = 1, and there are points at y = and ±Ay.) In the x direction, the grid 
starts near the axis at x min and extends out to some finite positive x max . We consider the 
two cases x min = — (v — \)Ax (staggered) and x m i n = —vAx (non-staggered). 

One can compute finite differences with a centered Cartesian 3D molecule in the interior 
plane of the slab, i.e. for all grid points at x > and y — 0, once boundary values are 
specified. The boundary values in the z direction and the positive x direction are assumed 
to be given as part of the general problem. As outlined above, the key point in the case of 
axisymmetry is that field values at points with y ^ 0, and also at the points with x < 0, can 
be obtained through a rotation and interpolation from field values on the half plane x > 
and y = 0. 

A. Rotation of Tensors 

Let us first discuss the rotation of the tensors and non-tensors that arise in typical numer- 
ical relativity computations. The basic formula, Eq. (|7|) below, is (aside from a few subtleties 



4 



y 



o/ o 
-6 — 6- 



□ 

-o 



ooooooon 



FIG. 1. A typical numerical grid in a constant z plane with Ax = Ay, no staggering, for a 
3D molecule of radius 1. The arcs show the rotation allowed by axisymmetry to obtain various 
field values on the slab boundary (circles). High order interpolation in the x direction may require 
several points at x < 0; for these points the rotation is shown by dashed arcs. The outer boundary 
at x max is given as part of the problem (squares). In the interior of the grid at < x < x max , 
y = (solid dots), the standard 3D finite difference stencil can be used. 



in its derivation) precisely what one expects for a vector in R 3 in Cartesian coordinates. Con- 
sider the rotation of a tensor field by an angle — 0o about the z axis. Equivalently, we can 
view this as keeping the field fixed, but rotating the coordinates in the opposite direction, 
i.e. rotating them by +0o- Taking this latter viewpoint, in cylindrical coordinates (p, 0, z) 
that are adapted to the rotation about the z axis, the rotated coordinates are 

p' = P, 0' = + 0o, z' = z, (3) 

or in the Cartesian coordinates given by p = \Jx 2 + y 2 , = arctan(y/x), 

x' = x cos 0o — 2/ sin 0(b y' — ^ sin O + ycos0 o , z' — z. (4) 

This coordinate transformation gives rise to the linear map 

fg x 'i\ ( cos^o -sin0 o \ 
(i2(0o)V> = (^j = I sin 0o cos 0o ol. (5) 

Note that i?(0 o )" 1 = R{-(j>o). 

Now consider arbitrary tensor fields T on a three-dimensional manifold E. Without 
referring to coordinates, a rotation (if it exists) defines a diffeomorphism R : S — > S, 
mapping points to points, and as such also defines a new tensor field R*T of the same 
contra- and covariant type as T. The diffeomorphism R is a symmetry transformation for 
the tensor field T if R*T = T. 

We assume, as is often the case in numerical relativity, that the domain of interest is 
covered by a single coordinate chart. The matrix of components of R* in the coordinate 
bases of the coordinate system (x 4 ) = (x,y,z) at a point p and the coordinate system 
(x n ) = (x', y', z') at the point R(p) equals the Jacobian matrix of the map R between the 
coordinates, which is given in (||), if we assume coordinates adapted to the rotation. The 
components of a tensor T at p, for example T iun '- 'j lt j 2 ,...(p)i transform according to 
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(R*Ty^>- jitj2 jR(p)) = /r, : /r, :! . . . (ir^yjr 1 )^ . . .T fc ^-- W2i ...(p). (6) 

In the case of axisymmetry (R*T = T) we therefore arrive at 

-,,,,...(.'••//. :) = R ll kl R i2 k 2 ■ ■ ■ {R- X ) h n {R~ l ) h J2 . . .T kl,k2, ~i lt i 2> _(\Jx 2 + y 2 , 0, z), (7) 

where R is given by @ with cos 4>o = x/p, sin O = y/p, and p = \/x 2 + y 2 . For example, for 
a vector T\ T l (x,y,z) = R l jT 3 (p,0, z). Equation (|7|) describes how to compute the values 
of a tensor at points outside the half-plane x > and y ^ from points within this plane. 

One type of non-tensor that is often used in numerical relativity is the partial derivative 
of tensors, 9fcT* 1 ' l2r " J1 j 2r ... Under coordinate transformations, the index k transforms by the 
chain rule of partial derivatives with the same Jacobian factor as appears in the coordinate 
transformation of a tensor, @. Using (||) and the product rule, there appear additional 
terms containing dj.R l j. However, using coordinates that are adapted to the axisymme- 
try (|) implies dkR l j = (in other words, R describes a rigid rotation), and hence the basic 
transformation law (0) is also valid for partial derivatives of tensors in this context. 

B. ID interpolation in the x direction 

We now turn to a brief discussion of the necessary interpolations. As made explicit 
in (|7p, we require fields at points (y/x 2 + y 2 , 0, z). As such points are not necessarily part of 
the numerical grid, a one-dimensional interpolation may be required. In this work we test 
Lagrange polynomial interpolation (e.g. [p7| , pS|1 ) with interpolation polynomials of degrees 2 
through 5, obtaining good results already with second order polynomials. Our evolution 
code uses second order finite difference molecules of size v = w = 1, but for interpolation 
polynomials of degree larger than 2 we need v > 1. Note that even the points at x = x max , 
can be obtained through interpolation by first applying the physical boundary for 
y = 0, then rotating to y ^ for x < x max , and then applying the physical boundary at 

This sort of interpolation is used in many areas of computational science (e.g. in adaptive 
mesh refinement), and thus should not be expected to pose a serious threat to the success of 
this technique. The principal issue of concern is whether interpolation might destabilize an 
evolution scheme that often is carefully crafted using particular forms of centered stencils 
and averages. This must be addressed on a case-by-case basis, either by Von Neumann or 
other analytical stability analysis, and/or by numerical experiment. 

Although we do not do it in practice, instead of explicitly computing and storing field 
values at the boundaries in the y direction, one can make the polynomial interpolation a 
part of the stencil, thereby explicitly reducing the 3D stencil to an inhomogeneous and 
asymmetric - but nonsingular - 2D stencil. In contrast, storing field values for a 3D stencil 
allows us to simply use numerical routines from existing Cartesian 3D codes. 

III. APPLICATIONS 

In this section we provide two important tests, both with the complete set of 3D, nonlin- 
ear Einstein evolution equations, performed with "Cartoon" in the Cactus code for numerical 
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relativity [j5|-0. However, we stress that the techniques developed here are applicable to 
many families of partial differential equations. That our tests are so successful, with such a 
complicated set of tensor PDEs as the 3 + 1 Einstein equations involving dozens of coupled 
nonlinear evolution equations, stresses this point. Simpler sets of equations can be expected 
to have fewer complications: less severe regularity issues at symmetry points (origin or axis), 
simpler transformation laws in the rotation to provide boundary conditions, etc. 



A. Einstein's Equations 



In this section we give a brief introduction to the basic equations we will be solving. The 
Einstein equations of general relativity in 3+1 form (see |||| for recent reviews and further 
references) are a complicated set of coupled, nonlinear partial differential equations for the 
symmetric tensor fields 7^ and (indices range from 1 through 3, so there are 12 field 
variables in all). The metric 7^ is the spatial part of the spacetime metric (which gives the 
invariant distance between two infinitesimally separated events): 

ds 2 = -(a 2 - P i Pi) dt 2 + 2pi dx { dt + 7ij dx { dx j . (8) 

The extrinsic curvature, K^, specifies how the t = constant slices are embedded in spacetime. 
The equations also involve the auxiliary tensor fields a and f3 l . These so-called "gauge" 
fields carry no dynamical information: they may be freely chosen for convenience. The lapse 
function a determines the proper time dr = adt measured by an observer falling normal 
to the time slice defined by t — constant. The shift vector (3 l determines the coordinate 
distance a constant-coordinate point moves away from the normal vector to the slice as one 
advances from one slice to the next. In the tests performed in this paper, we choose j3 l = 
for simplicity. 

The evolution equations can be written as: 

d tlij = -2aK {j + DiPj + DjPi, (9a) 



d t K i3 = -DiD ja + a[R i3 + {tr K)^ - 2K ik K k j \ 

+(3 k D k K tj + K lk Dj(3 k + K ]k DiP k . (9b) 

In these equations Rij is the 3-Ricci tensor, R is the 3-scalar curvature (both nonlinear 
functions of the metric 7^ and its first and second spatial derivatives), tri^ is the trace 
of Kij, and Di is the covariant derivative associated with the 3- metric 7^. With suitable 
choices for a and /3\ equations (|9|) are hyperbolic in 7^ and K^. 

The fields 7^ and are not completely freely specifiable: on each t = constant slice 
they must satisfy the four constraint equations 

H = R + (tr Kf - KijK ij = 0, (10a) 



H l = Dj(K ij - g ij tr K) = 0, (10b) 

where for later use we define the left hand side functions as H (known as the Hamiltonian 
constraint) and H\ Eqs. (0) are elliptic equations in 7^ and K^; in general they must be 
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solved numerically in order to obtain valid initial data ( f2~9f). However, one can show that 
once the constraints are satisfied on an initial slice, they are preserved by the evolution equa- 
tion (^|) (i.e. they stay satisfied for all future times). This statement holds for the continuum 
equations @ and fllOl) ; a finite difference evolution will in general only approximately satisfy 
the constraints at later times ( |3(|), and in fact the deviations of H and H % from zero are 
useful diagnostics of the evolution's numerical accuracy. 

To understand the physical meaning of these equations it is useful to consider an analogy 
with Maxwell's equations in a vacuum: 7^ and are analogous to the (vector) electric and 
magnetic fields E and B, the evolution equations ([|) are analogous to the Maxwell equations 
d t E = VxB and d t B = — V x E, and the constraint equations (10 ) are analogous to the 
Maxwell equations V-.E7 = V-B = 0, with H and H l analogous to V-J3 and V B. However, 
unlike Maxwell's equations, the Einstein equations are nonlinear, quite complicated (they 
have on the order of 1000 terms when written out in scalar form in terms of coordinate partial 
derivatives), and second order in space (though still first order in time). These properties 
make the Einstein equations difficult to treat numerically. 

As the main point of our paper is a technique for solving any set of evolution equations, 
we simply state that the particular form of the equations we use for this paper is a variant 
and ( |T0| ) recently put forth by Baumgarte and Shapiro JJT 



of 



of Shibata and Nakamura |52 



based on previous work 
We will refer to this as the BSSN formulation. It is to be 
noted that this formulation is quite similar in many respects to the Bona-Masso formula- 
tion [|33| -|35|| . In another paper [3f| we detail experiments carried out with these formulations 
on various spacetimes. As these details are not important to the results presented here, in 
this paper we focus only on the technique and its application to a very general class of partial 
differential equations, as represented by the Einstein equations. 



B. Black hole spacetimes 

In this section we report on the application of the Cartoon technique to the nonlinear 
evolution of black holes in the strong field regime. To illustrate the power of this technique, 
we focus on the case of a Schwarzschild black hole evolved with geodesic slicing. Geodesically 
sliced black hole evolutions have been used extensively to test black hole codes in ID |37| and 
2D in polar-spherical type coordinates, and in 3D Cartesian coordinates |T9| , |2T| , [22| , |24 . 



This black hole case has the advantages of (i) having an analytic and also a highly accurate 
numerical ID solution; (ii) providing a demanding test, as the slicing rapidly approaches the 
spacetime singularity, and hence spacetime curvature and metric functions grow rapidly and 
without bound until the code crashes at time t = irM, where M is the mass of the black 
hole; (iii) being one of the simpler testbeds in strong field regimes, as the geodesic slicing 
condition (a = 1) does not involve coupling additional elliptic or parabolic equations, which 
could complicate the demonstration of the technique. In the following section we consider a 
more complex system as a further test of the Cartoon method. 
The initial 3-metric for the Schwarzschild black hole is given by 

ds 2 = i>\dr 2 + r 2 (d6 2 + sin 2 6d<J) 2 )), (11) 

where the conformal factor is ip — (1 + jf). Here r is the isotropic radius, related to the 
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standard Schwarzschild radius r s by r s = (l + |^) 2 r. Transforming to Cartesian coordinates, 
we have 

ds 2 = ip 4 (dx 2 + dy 2 + dz 2 ), (12) 

where the Cartesian coordinates x, y, and z are related to the isotropic radius r in the usual 
way. The extrinsic curvature for the time symmetric, t — slice of this spacetime vanishes. 

We evolve this black hole in Cartesian coordinates, using the Cartoon technique. While 
typical 3D simulations of the Einstein equations are currently limited to grid sizes of roughly 
100 3 (ignoring the possibility of mesh refinements J38|), and the largest production simu- 
lations are currently at around 300 3 - 400 3 |Z5],|33| on a 128Gbyte SGI/Cray Origin 2000 
supercomputer with 256 processors, in this simulation we compute the evolution on a grid 
of N x x N y x N z = 1025 x 3 x 2049. (Our current implementation of the technique requires 
the z direction to range from +z max to — z max .) This is equivalent to a full 3D calculation 
of 2049 3 , over two orders of magnitude larger than the largest production simulations to 
date. Yet with this technique, the calculation requires only 0.07% of the memory required 
to do the full 3D simulation, and was run on 8 processors of an Origin 2000 over a period 
of 40 hours. 

In more detail, the present calculation was performed for a unit mass black hole, with 
Ax = Ay = Az = 0.01, utilizing an iterative-Crank-Nicholson (ICN) method of lines 
evolution scheme with three iterations, a time step At = 0.1, spatial gridresolution Ax = 
0.001, and 4th order Lagrange polynomial interpolation for the Cartoon algorithm. The total 
number of time steps for such a simulation was 3000 (9000 ICN iterations), corresponding 
to a time of t = 3.0. 

In Fig. |2| we show the metric function g xx on the x axis at selected times during the 
evolution. There is no hint of any instability in the result. At t — 3.0, a rather small 
difference between the numeric and analytic result (cmp. [F2l|,|T^]) can be seen at the peak of 
g xx . At t = 3.0 we also show g zz on the z axis, which at this scale is indistinguishable from 
the numerical g xx on the Ob 3jX1S ; clS it should be for this spherically symmetric data. The 
z direction is computed quite differently from the x direction, and suffers from the maximal 
inhomogeneity of the Cartoon stencil due to the interpolation near the z axis. Yet results 
along the 2-axis remain accurate and stable until the end of the simulation, giving a strong 
indication of the robustness of the technique. 

In Fig. [3] we demonstrate second order convergence of the Cartoon evolution. We use the 
Hamiltonian constraint if as a diagnostic of the evolution's numerical accuracy: since H 
vanishes analytically, it should converge to zero as the grid resolution increases. In particular, 
for second order finite differencing, each time the grid resolution is doubled H should shrink 
by a factor of 4 ( |3^]). The figure shows the maximum value of H in the grid as a function 
of time, plotted for the 3 different grid resolutions Ax = Ay = Az = 0.02, 0.04, 0.08, with 
each successive plot divided by 1,4, 16 respectively, so the values should be identical at the 
different resolutions. Second order convergence is directly evident from the figure. 

It is to be noted that not only is such a simulation impossible today in full 3D due 
to memory constraints, but the 2D Cartoon run has far higher resolution than has been 
achieved even in 2D codes designed to treat similar problems. Axisymmetric simulations of 
black holes, usually performed in spherical-polar coordinates, are typically not performed 



with more than 300 radial and 50 angular zones [40,|]]mE2|. Inherent axis instabilities 
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FIG. 2. We show results for a geodesically sliced Schwarzschild black hole, evolved in 
3D with the Cartoon technique. The metric function g xx is shown on the x axis at times 
t/M = 0.0, 0.5, . . . , 3.0, where M is the mass of the black hole. At t/M = 3.0, g zz on the z axis is 
also plotted, but it is indistinguishable from g xx as expected for spherical symmetry. Both agree 
well with the analytic result. 
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FIG. 3. The maximum of the Hamiltonian constraint for a geodesically sliced black hole at 
three different resolutions. Each graph is scaled proportionally to the inverse square of the grid 
resolution. The equality of the rescaled values for the different resolutions indicates second order 
convergence. 



11 



usually create difficulties in simulations with higher resolution than this, in part because 
angular zones closer to the axis have more delicate regularity behavior to enforce. It is 
precisely these sorts of problems Cartoon is designed to overcome. 

The problem with standard axisymmetric codes are exacerbated unless certain gauges are 
used that force troublesome terms in the three-metric to vanish. For example, as discussed 
in Ref. || , a vicious axis instability arose in simulations of spherical or distorted black holes 
in axisymmetry, causing a rapid code crash, until a gauge condition was developed that 
diagonalized the three-metric. This condition required the solution of an elliptic equation 
to determine j3 l at each time step, which is very time consuming to solve numerically. 
Furthermore, even with this gauge choice, if the code is run with too high a resolution, 
axis instabilities are again encountered, and these instabilities eventually crash the code at 
late times. Similar problems with axis instabilities have been encountered in the case of 
rotating and colliding black holes @,g|5],|6|,0,||,|3 . 

As reported recently in Ref. ||24j| , a case similar to the one presented here, a geodesically 
sliced perturbed black hole, was studied to compare results from a 3D code in Cartesian 
coordinates to a traditional 2D axisymmetric code. In order to compare metric functions 
directly, the codes had to be run with the same spatial gauge, which had vanishing shift. 
Although the agreement was excellent for as long as it could be computed, due to the axis 
instability the axisymmetric code crashed far earlier than the time at which the slice hit 
the singularity. On the other hand, the 3D Cartesian code, due to its lack of coordinate 
singularities, was able to run accurately all the way to the singularity. With our Cartoon 
technique, we are able to do the same simulation in Cartesian coordinates, with the same 
accuracy and stability, but with far less memory than is required in full 3D, and with far 
more stability than is possible in a standard axisymmetric code. 



C. Gravitational wave spacetimes 

While the above simulation of time slices approaching a black hole singularity provide a 
strong test of the techniques we have developed, the underlying system is spherically sym- 
metric, and does not contain gravitational waves. In order to further test and confirm this 
technique, we now turn to a completely different spacetime system, one not initially contain- 
ing no black hole. Instead, the system we choose is that of highly nonlinear gravitational 
waves. 

Gravitational waves are often considered in the linearized regime, where they are small 
disturbances on some background (often flat) spacetime that propagate at the speed of 
light. However, as the Einstein equations are nonlinear, for cases where the waves have 
sufficiently large amplitude they can affect the background spacetime on which they prop- 
agate. Even more, for strong enough waves there is no background spacetime; the waves 
have to be treated fully nonlinearly, and under extreme conditions they can even collapse in 
on themselves, forming a black hole when none existed previously. 

Low amplitude gravitational waves have provided testbeds of 3D numerical relativity 
codes for the last decade, and even there they have provided a strong challenge H32| , ^0| . [l7| . 
Extreme gravitational wave simulations, where the waves form their own background, and 
where they may actually form black holes, have only been possible to evolve in full 3D 
codes during the last year HSfl. For waves above a certain critical amplitude a black hole 
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forms, while for waves just below critical a rich pattern of oscillations develops as the waves 
teeter on the edge of forming a black hole, and then eventually they disperse and the system 
returns to flat space (in a highly nontrivial coordinate system). 

In this paper we use one such simulation to illustrate the strength of the Cartoon tech- 
nique, and compare with a full 3D simulation as a measure of its accuracy. As our main 
purpose here is to demonstrate that the technique is robust, even under very demanding 
simulations, we will not go into detail of the physics of these simulations. The interested 



reader is asked to consult Ref. [48| for more details 



As in Ref. 



IS |, we take as initial data a pure gravitational wave data set, based on the 
axisymmetric ansatz of Brill |EJ , and later studied by Eppley P5|J50| , |5T| and others |52| , |53 
The metric takes the form 



ds 2 



e 2 « (dp 2 + dz 2 ) + 



P 



ty 4 ds 



(13) 



where q is a free function subject to certain boundary conditions. We consider a function q 
of the form 



q 



a p 2 e 



(14) 



where a is a constants, and r 2 



p 2 + z 2 



sec 



, x X Ny 



The same system is solved 
x N g = (129 x 129 x 129) 



HI for 2D, and @g| for full 3D Cartesian 
coordinates). We consider data with the amplitude a = 4 which corresponds to a strong, 
axisymmetric, equatorial plane symmetric gravitational wave. We choose the extrinsic cur- 
vature to vanish (time symmetric initial data). As shown in Ref. [[ISj], this initial data 
collapses in on itself initially, but after a series of reverberations it disperses, leaving flat 
space in its wake. 

Taking this form for q, we solve the Hamiltonian constraint equation (|10a| ) numerically 
using a multigrid elliptic solver on an N x x N y x N z = (129 x 3 x 257) grid, with Ax = Ay = 
Az = 0.04, and the outer boundary at 5.12. As this axisymmetric initial data set is also 
symmetric about all coordinate planes, it is possible to evolve it as a full 3D system in just 
one octant (x,y,z non- negative), as was also done in Ref. 
numerically in full 3D (with octant symmetry) on an N, 
grid, again with Ax = Ay = Az = 0.04, as a comparison simulation. 

Both the full 3D system and the slab system were evolved to a time when the gravitational 
waves had largely left the system (through outgoing boundary conditions applied on the 
evolved function; see Ref. |18j for details). In Fig. f|, as a sensitive measure of the evolution 
we show the minimum value of the lapse function a as a function of time throughout the 
evolution. We have found this to be a good indicator of the evolution of the system pSfl . 
and also, as it sits on the axis and origin of the system, it should be especially sensitive 
to any problems that arise during the evolution due to the Cartoon procedure. We show 
results for 3D, and for Cartoon using interpolation of order 2, 3, and 4. The agreement 
is excellent. There is a discernible difference between the 3D and Cartoon runs after time 
t = 10, starting with a slight bump that can be seen in the slab evolution that is not present 
in the full 3D simulation. This is due to slight differences in the boundary treatment in 
the two cases. There is no known perfect outgoing wave condition for general relativity 
(nor even a clear way to identify what a wave is in a nonlinear situation such as this), and 
small differences in boundary treatment can lead to differing amounts of reflection at the 
boundaries. The different geometries of the two simulations lead to different characteristics 
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FIG. 4. The minimum of the lapse during the evolution of a strong Brill wave, which generically 
occurs at the origin and is a sensitive measure of the spacetime evolution. We compare three 
Cartoon runs using interpolation order 2, 3, and 4, with a 3D run. 

of the boundary treatment, which are ultimately reflected in the results. The origin of this 
difference is well-understood, to be expected, and minor, and is unrelated to the main thrust 
of this paper. 



We have developed a novel technique for computational simulations with certain sym- 
metries. This Cartoon technique allows simulations to be performed in 3D Cartesian coor- 
dinates, thereby avoiding difficulties associated with coordinate singularities that often lead 
to numerical instabilities. Compared to the full 3D Cartesian approach, Cartoon allows a 
huge savings in both memory and computational time, without introducing problems com- 
monly associated with the singular coordinate systems which are usually employed. The 
Cartoon technique should be useful for any system of partial differential equations, and we 
have shown that it works well in one of the most complicated sets of equations in theoretical 
physics, Einstein's equations of general relativity. 

While this approach to axisymmetric simulation inherits many of the advantages of 
3D Cartesian simulations, there are also some disadvantages which carry over. For example, 
a 2D code using angular coordinates can align spherical boundaries at constant coordinate 
value. Furthermore, a logarithmic or other stretched radial coordinate can be introduced 
that significantly improves resolution for problems with asymptotic 1/r fall-off. Even though 
neither of these techniques are readily incorporated in the Cartoon method, we nonetheless 



IV. CONCLUSIONS 
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expect Cartoon to have a great impact on numerical relativity, as general purpose stable 2D 
codes are not currently available in the field. 

Furthermore, as most fully 3D simulations are carried out in the same Cartesian coordi- 
nates as are found in Cartoon, experience gained from the accelerated Cartoon simulations 
should carry over directly to the full 3D work. This has generally not been the case previ- 
ously, as special coordinate systems used for axisymmetric or spherically symmetric systems 
also required special treatments and gauge conditions that simply were not applicable in 
3D Cartesian coordinates. Hence with our new technique systems can be studied with the 
same coordinate systems, with the same gauges, and with the same analysis tools as they 
will be when performed in full 3D. Finally, this technique has the potential to allow for 
a number of intrinsically axisymmetric systems to be studied with unprecedented stability 
and accuracy, since simulations with resolutions of several thousand grid points on a side 
are readily achievable. 

We have implemented this technique in the Cactus code for 3D numerical relativity and 
the Cactus Computational Toolkit, and we expect that it will be a powerful addition to the 
field of numerical relativity. 
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